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Aegean  Seas”  by  Kathryn  A.  Kelly 

A  column  mixed  layer  model  is  run  for  the  Aegean  and  eastern  Mediterranean  Seas  and  results  are 
compared  with  SST  and  for  the  autumns  of  2001-2004  to  determine  whether  errors  in  modeling  SST 
can  be  attributed  to  missing  model  physics,  primarily  mixing.  Discrepancies  between  the  modeled  and 
observed  SST  were  found  to  be  about  the  same  size  as  differences  between  SST  products,  so  that 
model  errors  could  not  be  definitively  isolated.  However,  much  canbe  learned  about  the  accuracy  of 
the  forcing  and  observed  fields  from  this  study.  SST  fluctuations  with  temporal  scales  of  10—14  days 
were  determined  to  be  caused  by  latent  heat  flux  variability,  which  is  somewhat  underestimated  by  the 
second  NCEP/NCAR  Reanalysis  (NCEP2)fluxes.  Spatial  variations  in  air-sea  fluxes  on  spatial  scales 
not  resolved  by  the  Reanalysis,  as  determined  by  comparisons  between  NCEP2  and  scatterometer 
winds,  substantially  degrade  model  performance.  Contributions  from  advection  are  relatively  small. 
Biases  in  the  NCEP2  precipitation-minus-evaporation  fields  cause  systematic  errors  in  mixed  layer 
depth  and  in  temperature.  Relaxing  the  salinity  to  climatological  values,  combined  with  small  bias 
corrections  to  air-sea  fluxes  gives  better  agreement  with  observed  SST.  Model  performance  was 
judged  by  comparing  observed  and  modeled  temperature  tendency,  a  more  stringent  comparison  than 
with  SST;  good  results  were  obtained  in  the  eastern  Mediterranean  and  Aegean  Seas,  but  not  in  the 
central  Mediterranean.  A  relatively  poor  match  with  an  infrared-based  SST  product  is  apparently 
caused  by  cloud-contamination,  not  the  mixed  layer  model,  a  result  established  by  comparisons  with 
microwave  SST  and  cloud  cover. 
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Abstract.  A  column  mixed  layer  model  is  run  for  the  Aegean  and  eastern 
Mediterranean  Seas  and  results  are  compared  with  SST  and  for  the  autumns  of 
2001-2004  to  determine  whether  errors  in  modeling  SST  can  be  attributed  to 
missing  model  physics,  primarily  mixing.  Discrepancies  between  the  modeled 
and  observed  SST  were  found  to  be  about  the  same  size  as  differences 
between  SST  products,  so  that  model  errors  could  not  be  definitively  isolated. 
However,  much  can  be  learned  about  the  accuracy  of  the  forcing  and  observed 
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somewhat  underestimated  by  the  second  NCEP/NCAR  Reanalysis  (NCEP2) 
fluxes.  Spatial  variations  in  air-sea  fluxes  on  spatial  scales  not  resolved 
by  the  Reanalysis,  as  determined  by  comparisons  between  NCEP2  and 
scatterometer  winds,  substantially  degrade  model  performance.  Contributions 
from  advection  are  relatively  small.  Biases  in  the  NCEP2  precipitation- 
minus-evaporation  fields  cause  systematic  errors  in  mixed  layer  depth  and  in 
temperature.  Relaxing  the  salinity  to  climatological  values,  combined  with 
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SST.  Model  performance  was  judged  by  comparing  observed  and  modeled 
temperature  tendency,  a  more  stringent  comparison  than  with  SST;  good 
results  were  obtained  in  the  eastern  Mediterranean  and  Aegean  Seas,  but  not 
in  the  central  Mediterranean.  A  relatively  poor  match  with  an  infrared-based 
SST  product  is  apparently  caused  by  cloud-contamination,  not  the  mixed 
layer  model,  a  result  established  by  comparisons  with  microwave  SST  and 
cloud  cover. 
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1.  Introduction 

Sea  surface  temperature  (SST)  variations  are  a  com¬ 
plicated  response  to  air-sea  heat  and  freshwater  fluxes, 
wind,  mixing,  and  advection.  In  a  region  as  complex 
geographically  as  the  Mediterranean  Sea  one  expects 
these  processes  to  have  different  balances  in  different 
regions.  In  particular,  higher  levels  of  mixing  associ¬ 
ated  with  internal  waves,  tides,  or  other  geographically 
correlated  processes  may  produce  temperature  anoma¬ 
lies  that  are  not  predicted  by  simple  mixed  layer  mod¬ 
els.  Conversely,  systematic  errors  between  observed  and 


modeled  temperatures  may  indicate  the  location  and 
extent  of  anomalous  mixing  processes.  More  generally, 
given  sufficiently  accurate  forcing  fields  and  observa¬ 
tions  for  comparison,  the  question  addressed  here  is: 
Do  temperature  errors  indicate  missing  model  physics? 

This  study  was  motivated  by  a  workshop  on  the 
Aegean  Sea,  sponsored  by  the  U.S.  Office  of  Naval  Re¬ 
search  [Sofianos  et  al,  2002],  and  a  field  program,  of 
which  one  of  the  goals  is  to  understand  the  contribu¬ 
tions  of  tides  and  internal  waves  to  ocean  mixing.  Par¬ 
ticularly  in  the  fall  when  thermal  stratification  is  large, 
internal  waves  and  tides  would  be  expected  to  increase 
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vertical  mixing  and  to  decrease  SST.  Related  scientific 
questions  raised  in  the  workshop  report  include:  How 
is  the  wind  stress  pattern  affected  by  the  presence  of 
complex  orography?  How  do  the  resulting  patterns  of 
surface  wind  stress  around  the  islands  affect  eddies,  in¬ 
ternal  waves,  and  mixing?  What  are  the  patterns  and 
dominant  processes  affecting  oceanic  convection  and  re- 
stratification?  How  are  these  processes  modulated  by 
islands,  plateaux,  and  eddies?  While  these  questions 
can  only  be  conclusively  answered  in  conjunction  with 
extensive  field  measurements  and  concurrent  modeling 
studies,  an  analysis  of  the  increasingly  accurate  forcing 
fields  and  observations  currently  available  may  give  a 
useful  overview  and  assist  in  designing  more  detailed 
studies. 

Several  large-scale  processes  that  likely  influence  tem¬ 
peratures  in  different  regions  of  the  Mediterranean  and 
Aegean  Seas  have  been  investigated  using  field  pro¬ 
grams  and  numerical  studies.  Sea  level  has  been  rising 


servations  available  for  comparison  in  the  Aegean  Sea, 
as  well  as  some  visible  images  of  the  region  from  the 
MODIS.  Analyses  of  observations  showing  regions  of 
high  variability  are  presented  in  Section  2.  Section  3 
gives  the  methodology,  including  the  model  set  up,  and 
the  adjustments  to  the  forcing  fields.  In  Section  4  the 
model  and  observations  are  compared,  and  the  sensitiv¬ 
ity  of  the  model  to  forcing  is  examined.  The  summary 
and  conclusions  are  presented  in  Section  5. 

2.  Evaluation  of  Observations  and 
Forcing 

Several,  primarily  satellite,  data  sets  were  examined 
to  provide  information  on  regional  variability,  including 
the  existence  of  internal  waves,  the  effects  of  topogra¬ 
phy  on  winds,  and  temperature  variations.  Both  tem¬ 
perature  and  wind  accuracies  were  evaluated  by  com- 
cuiia  cutu  iiumciicai  atuuico.  oca  icvci  uaa  uccu  ik>ui£  parisons  with  observations  from  a  network  of  buoys  in 
in  the  Mediterranean  Sea,  as  observed  by  the  TOPEX/Posei^RR  Aegean  Sea,  maintained  by  the  Hellenic  Centre  for 
radar  altimeter  [Cazenave  et  al.,  2001],  particularly  in  Marine  Research, 
the  region  east  of  Crete;  SST  and  hydrographic  data 
show  corresponding  increases  in  water  temperature. 

A  recent  study  by  Fukumori  [2006]  showed  that  the 
Mediterranean  Sea  has  a  relatively  rapid  and  nearly 
uniform  barotropic  (seiche-like)  response  to  wind  stress 
with  time  scales  of  weeks  to  months,  based  on  model 
experiments  and  altimeter  data.  A  possible  influence 
of  the  wintertime  North  Atlantic  Oscillation  has  been 
found  in  the  circulation  of  the  Western  Mediterranean 
Sea  [Vignudelli  et  al .,  1999],  based  on  a  time  series  of 
currents  in  the  Corsica  Channel.  A  numerical  simu¬ 
lation  of  circulation  by  Pinardi  et  al.  [1997]  showed 
that  the  Eastern  Mediterranean  Sea  has  considerably 
more  interannual  variability  than  the  Western  Mediter¬ 
ranean.  Unlike  in  the  west,  strong  wintertime  winds 
and  heat  fluxes  in  the  Eastern  Mediterranean  can  mod¬ 
ify  the  ocean  circulation  and  structure  to  overcome  the 
normal  seasonal  cycle  the  following  summer.  Quasi¬ 
stationary  and  recurrent  eddies  have  been  observed  in 
the  eastern  Mediterranean  Sea  using  aerial  surveys  and 
drifting  buoys  [Matteoda  and  Glenn ,  1996]. 

Against  this  background  of  regional  variations  in 
ocean  and  forcing  variability,  the  ability  of  a  well-known 
mixed  layer  model  to  simulate  sea  surface  temperatures 
(SST)  is  evaluated  for  the  eastern  Mediterranean  Sea 
(Levantine  and  Ionian  Basins)  and  for  the  Aegean  Sea. 

The  Price- Weller-Pinkel  [Price  et  a/.,  1986]  mixed  layer 
model  is  used  here,  forced  by  new  air-sea  flux  and  wind 
products.  Several  air-sea  flux  products  are  evaluated 
to  determine  their  effectiveness;  the  model  results  are 
then  compared  with  two  new  SST  products.  In  addi¬ 
tion,  there  are  some  in  situ  temperature  and  wind  ob- 


2.1.  Satellite  Observations  of  Internal  Waves 

To  determine  the  effect  of  internal  waves  on  ocean 
mixing,  it  would  be  useful  to  have  a  statistical  descrip¬ 
tion  of  when  and  where  internal  waves  occur.  Inter¬ 
nal  waves  have  been  observed  using  a  variety  of  sensors 
[Global  Ocean  Associates ,  2004],  including  an  imaging 
radar  on  the  Space  Shuttle.  Two  sources  of  data  exist 
from  which  a  statistical  description  might  be  obtained: 
high-resolution  visible  images  and  synthetic  aperture 
radar  (SAR)  images. 

True  color  (visible)  images  from  MODIS  (Figure  la) 
on  either  the  Aqua  or  Terra  satellites  are  readily  avail¬ 
able  with  a  resolution  of  about  250m  from  NASA’s  God¬ 
dard  Space  Flight  Center.  Modulations  of  the  ocean 
surface  roughness  can  only  be  seen  when  highlighted  by 
sunglint  in  cloud-free  regions;  sun-ocean  surface  view¬ 
ing  geometry  limit  the  useful  observation  period  to  the 
months  of  June,  July,  and  August.  All  available  data 
(64  images)  from  both  platforms  for  2002-2004  were  ex¬ 
amined.  During  this  period,  there  were  no  obvious  ex¬ 
amples  of  internal  waves  in  the  eastern  Mediterranean 
or  Aegean  Seas;  however,  there  are  no  objective  ways  to 
determine  whether  the  observed  signatures  are  caused 
by  the  atmosphere  or  by  processes  internal  to  the  ocean. 

Browse  images  of  SAR  data  from  the  European  Space 
Agency  were  also  examined.  SAR  data  is  not  restricted 
by  viewing  geometry  or  clouds.  Again,  there  were  no 
clear  examples  of  internal  waves  in  the  SAR  images,  and 
recent  studies  that  suggest  that  no  algorithm  exists  to 
extract  the  signature  of  internal  waves  in  the  presence  of 
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many  other  types  of  oceanographic  phenomena  [ Ivanov , 

2002]. 

Nevetheless,  the  MODIS  images  provide  useful  infor¬ 
mation  about  variability  in  the  region.  Bright  patches 
are  observed  (Figure  la)  on  the  south  side  of  Crete  and 
the  numerous  islands  in  the  Aegean  Sea  in  about  70% 
of  the  images, usually  during  periods  of  strong  (south¬ 
ward)  winds.  Wave  or  eddy-like  surface  patterns  occur 
in  about  50%  of  the  images,  usually  in  a  location  near 
the  edge  of  the  sunglint,  suggesting  that  these  features 
may  be  obscured  by  the  bright  patches  that  appear  un¬ 
der  high-wind  conditions. 


2.2.  Wind  Fields 

High-resolution  data  from  the  SeaWinds  scatterom- 
eter  on  QuikSCAT  (Figure  lb)  reveal  topographically 
modified  wind  fields  that  resemble  wakes  on  the  down¬ 
wind  side  of  islands  [Chelton  et  a/.,  2004].  The  bright 
patches  in  MODIS  images  (Figure  la)  likely  represent 
a  sea  surface  response  to  the  winds:  wind  shielding  by 
island  topography  results  in  a  smoother  ocean  surface 
that  in  turn  reflects  more  sunlight,  producing  the  bright 
patches  in  the  MODIS  images.  Such  high-resolution 
wind  fields  offer  the  potential  for  improved  modeling  of 
ocean  processes.  Correlations  between  the  buoy  winds 
and  collocated  12.5km-resolution  QuikSCAT  winds  are 
typically  0. 8-0.9.  To  force  an  ocean  model,  the  origi¬ 
nal  swath-oriented  wind  vectors  must  be  gridded  at  a 
lower  resolution,  consistent  with  the  revisit  time  of  the 
satellite  (about  12  hours).  At  the  time  of  this  study, 
the  high-resolution  QuikSCAT  winds  shown  here  were 
only  available  for  calendar  year  2003.  Therefore,  the 
winds  used  in  the  modeling  study  were  the  standard 
25-km  winds,  objectively  mapped  to  a  l°x  1°  grid  with 
approximate  4-day  resolution  (to  maintain  consistent 
data  quality,  Kelly  et  al.  [1999]).  Despite  the  high  de¬ 
gree  of  smoothing,  correlations  with  buoy  winds  do  not 
drop  dramatically  with  gridding:  correlations  between 
buoy  and  the  gridded  winds  are  typically  0.75-0.8.  A 
detailed  evaluation  of  the  high-resolution  winds  and  the 
effects  of  gridding,  relative  to  the  winds  from  the  buoys, 
will  be  the  subject  of  a  future  report. 

The  gridded  QuikSCAT  l°x  1°  wind  products  have 
considerably  higher  spatial  resolution  than  most  cur¬ 
rently  available  numerical  weather  prediction  products, 
such  as  NCEP2  (see  Section  4.3).  While  improved  res¬ 
olution  in  wind  stress  and  wind  stress  curl  would  likely 
give  improvements  in  dynamical  models,  wind  stress 
makes  a  relatively  small  contribution  to  the  variabil¬ 
ity  of  the  mixed  layer  model,  for  which  the  dominant 
forcing  is  the  air-sea  flux.  A  recent  study  by  Jiang 
et  al.  [2005]  suggests  that  combining  QuikSCAT  winds 


Figure  1.  Island  wakes  on  the  leeward  (south)  side 
of  Crete  on  28  July  2003.  (a)  MODIS  true-color  im¬ 
age  from  Aqua  at  1120  UTC  and  (b)  high- resolution 
QuikSCAT  wind  vectors  at  1700  UTC.  Color  contours 
are  wind  speed  with  contour  interval  of  1  m  s"1.  Low 
wind  speeds  in  (b)  correspond  to  lighter  colors  in  (a). 
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and  satellite  SST  with  lower-resolution  air  temperature 
and  humidity  produces  an  improved  flux  product  in  the 
equatorial  Pacific,  where  latent  heat  fluxes  clearly  dom¬ 
inate,  but  this  may  not  be  true  for  other  regions.  To 
take  full  advantage  of  the  improved  resolution  in  winds 
(and  SST)  for  heat  fluxes,  air  temperature  and  humid¬ 
ity  would  also  need  to  be  available  at  higher  spatial 
resolution.  An  estimate  of  the  magnitude  of  the  errors 
induced  by  poor  spatial  resolution  in  the  flux  fields  is 
included  in  Section  4.3,  using  the  region  south  of  Crete. 

2.3.  Sea  Surface  Temperature  (SST) 

Two  SST  products  are  used  in  the  evaluation  of 
the  model  performance  (Figure  2).  The  first  is  a  rel¬ 
atively  new  product  from  the  National  Centers  for  En¬ 
vironmental  Prediction  (NCEP),  designated  Real  Time 
Global  (RTG).  Although  it  is  derived  from  the  same  in¬ 
frared  sensors  that  are  used  for  the  standard  “Reynolds” 
[Reynolds  et  al,  1994]  SST,  the  spatial  resolution  of 
the  RTG  SST  is  much  higher.  The  second  SST  prod¬ 
uct  used  here  is  an  optimally  interpolated  version  of 
microwave  SST  (from  the  TRMM  Microwave  Imager) 
available  from  Remote  Sensing  Systems.  Because  mi¬ 
crowave  sensors  can  measure  through  clouds,  unlike  in¬ 
frared  sensors,  the  microwave  SST  products  axe  poten¬ 
tially  more  accurate,  particularly  in  regions  of  persistent 
cloud  cover.  Disadvantages  of  the  microwave  sensors 
are  that  it  has  inherently  lower  spatial  resolution  («  50 
km)  and  that  any  land  within  its  relatively  large  field- 
of-view  renders  the  data  unusable.  Thus,  microwave 
SST  is  not  reliable  in  most  of  the  Aegean  Sea  or  near 
any  coastline. 

Before  evaluating  the  ability  of  a  model  to  reproduce 
observed  mixed  layer  temperatures,  the  SST  products 
were  first  compared  with  temperatures  from  buoys  in 
the  Aegean  Sea  at  3  meters  below  the  surface  at  3-hr 
intervals  (See  Figure  3  for  buoy  locations).  The  com¬ 
parison  time  period  was  from  January  2001  to  mid-May 
2004.  Microwave  SST  was  available  at  only  the  two 
southernmost  buoys.  For  these  two  buoys,  RTG  SST 
had  a  negligible  bias  of  about  0.2°C,  whereas  TMI  SST 
had  a  bias  of  approximately  1.5°C.  Correlations  be¬ 
tween  SST  and  buoy  temperatures  were  examined  after 
first  removing  the  seasonal  cycle  (once  and  twice  per 
year  harmonics);  correlations  for  all  SST  series  were 
statistically  significant  with  RTG  values  of  0.50  and 
0.49  and  TMI  values  of  0.61  and  0.46.  The  proximity 
of  land  throughout  the  Aegean  Sea  likely  degrades  the 
TMI  SST  in  the  buoy  comparisons  and  causes  a  warm 
bias.  In  the  central  Mediterranean,  TMI  appears  to  be 
more  accurate  than  RTG  (see  Section  4.1,  below). 

The  quantity  used  to  determine  how  well  the  PWP 


(a)  TMI  SST  for  Nov  15.  2002 


(b)  RTG  SST  and  buoy  locations  (*) 
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Figure  2.  Maps  of  SST  products  for  17  October  2001. 
SST  maps  derived  from  (a)  the  microwave  sensor  TMI 
and  (b)  the  high-resolution  NCEP  RTG  product.  Loca¬ 
tions  of  three  buoys  in  the  Aegean  Sea  used  to  evaluate 
satellite  fields  shown  as  asterisks.  Buoys  are  named 
after  nearby  islands:  from  south  to  north,  Avgo,  San¬ 
torini,  and  Mykonos. 


model  performed  was  actually  temperature  tendency, 
dT/dt,  rather  than  temperature  itself  because  it  is  the 
rate  of  temperature  change  that  is  directly  related  to 
the  forcing.  Therefore,  this  quantity  was  also  compared 
between  SST  and  the  3-m  temperatures  (Figure  4).  A 
slight  modification  was  necessary  for  the  TMI  data  to 
increase  the  accuracy  of  dT / dt .  Because  the  mapped 
TMI  data  are  used  in  near  real  time,  the  interpolation 
algorithm  uses  only  data  from  previous  times.  In  the 
relatively  rare  event  of  missing  data  owing  to  rain,  SST 
from  the  previous  clear  period  is  used,  resulting  in  an 
unrealistically  high  occurrence  of  zeroes  in  the  value  of 
dT/dt.  Therefore,  wherever  dT/dt  was  found  to  be 
exactly  zero,  the  temporally  constant  SST  values  were 
replaced  by  a  linear  interpolation  between  previous  and 
subsequent  SST  estimates. 

To  estimate  a  sensible  value  for  temperature  ten¬ 
dency,  it  is  necessary  to  smooth  the  SST  fields  tempo¬ 
rally  to  reduce  errors  in  the  SST  products,  while  retain¬ 
ing  the  highest  frequency  variations  actually  resolved  by 
the  fields.  To  determine  an  appropriate  temporal  low- 
pass  filter  to  apply  to  the  SST  data,  correlations  were 
computed  between  temperature  tendency  from  SST  and 
from  daily  averaged  3-m  buoy  temperatures,  after  ap¬ 
plying  a  series  of  lowpass  filters  with  half-power  peri¬ 
ods  from  4  to  12  days.  Correlations  between  the  two 
filtered  series  of  dT /dt  increased  steeply  up  to  periods 
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(a)  35.6N,  25.6E 


(c)  37. 5N,  25.5E 
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Figure  3.  Comparison  of  SST  products  with  buoy  tem¬ 
peratures.  Time  series  at  the  three  southernmost  buoy 
locations  shown  in  Figure  2a  for  the  second  half  of  2002. 
TMI  SST  (thick  line),  RTG  SST  (thin  line)  and  3-m 
buoy  temperatures  (dots). 


(a)  unfiltered  dT/dt  at  36.3N.  25.5E 


(b)  filtered  dT/dt  at  36.3N,  25. 5E 


Figure  4.  Temperature  tendency  of  SST  and  buoys, 
(a)  Time  series  of  dT/dt  at  Santorini  (compare  Fig¬ 
ure  3b)  for  the  second  half  of  2002  for  RTG  SST  (solid 
line)  and  3-m  buoy  temperatures  (dots),  (b)  Lowpass 
filtered  tendency  for  RTG  SST  with  half-power  of  9 
days.  Tendency  from  buoys  repeated  for  comparison. 


of  approximately  8-10  days  (indicating  a  reduction  in 
SST  errors)  and  then  remained  relatively  flat  with  in¬ 
creasingly  longer  periods.  This  analysis  suggests  that 
SST  fluctuations  with  periods  less  than  about  9  days 
are  dominated  by  errors  in  the  SST  maps  or  are  asso¬ 
ciated  with  the  skin  temperature  of  the  ocean.  There 
are  energetic  fluctuations  in  the  buoy  temperatures  with 
periods  shorter  than  9  days,  but  these  fluctuations  are 
not  resolved  by  the  SST  products.  In  comparisons  with 
the  model,  all  SST  data  have  been  lowpass  filtered  with 
a  half-power  period  of  9  days. 

2.4.  Air-Sea  Heat  Fluxes 

The  model  was  forced  by  air-sea  heat  fluxes  and 
precipitation-minus-evaporation  (or  freshwater  fluxes), 
primarily  from  the  second  NCEP/NCAR  Reanalysis 
(NCEP2)  [Kistler  et  al .,  2001].  Different  combinations 
of  radiative  and  turbulent  heat  fluxes  were  used  to  force 
the  model,  including  both  NCEP2  radiative  and  turbu¬ 
lent  heat  fluxes,  and  radiative  fluxes  from  the  Inter¬ 
national  Satellite  Cloud  Climatology  Project  [Zhang  et 
al .,  2003].  An  additional  set  of  latent  and  sensible  heat 
flux  fields  were  derived  by  using  the  NCEP2  variables  in 
the  COARE  algorithm  [Fairall  et  al .,  1996],  version  3.0. 
The  PWP  model  also  requires  wind  stress  (for  horizon¬ 
tal  momentum)  and  wind  stress  curl  for  calculating  ver¬ 
tical  velocity  and  vertical  diffusion.  Gridded  QuikSCAT 
wind  stress  and  wind  stress  curl  fields  were  used  to  force 
the  model;  however,  as  discussed  in  Section  4.2,  the  re¬ 


sults  are  relatively  insensitive  to  the  accuracy  of  these 
products. 

Several  heat  flux  products  were  compared  and  tested 
for  their  ability  to  simulate  mixed  layer  temperature  in 
the  PWP  model.  The  flux  components  are  first  exam¬ 
ined  separately:  radiative  fluxes  (shortwave  and  long¬ 
wave)  and  turbulent  fluxes  (latent  and  sensible).  The 
oceanographic  convention  is  used  here,  that  is,  posi¬ 
tive  fluxes  represent  heat  fluxed  into  the  ocean.  Radia¬ 
tive  fluxes  include  daily  fields  from  NCEP2  and  daily 
fields  from  ISCCP.  Turbulent  fluxes  include  the  NCEP2 
product  as  well  as  a  flux  product  derived  by  using 
daily  NCEP2  atmospheric  variables  in  the  COARE  al¬ 
gorithm. 

Mean  values  over  2001-2004  were  compared  for  each 
flux  component  and  product.  Shortwave  radiative  fluxes 
from  NCEP2  are  on  average  about  20Wm-2  stronger 
than  the  ISCCP  fluxes  (Figure  5),  which  incorporate 
cloud  forcing  from  observations.  Longwave  fluxes  (not 
shown)  are  comparable  in  magnitude.  NCEP2  turbu¬ 
lent  fluxes  are  also  compared  with  the  NCEP2/COARE 
fluxes.  Latent  heat  flux  products  (Figure  6)  differ  by 
about  20Wm-2  overall,  with  NCEP2  more  negative, 
except  near  the  northern  coast  of  Africa,  where  the 
influence  of  land  appears  to  make  the  NCEP2  fluxes 
change  abruptly  to  values  less  negative  than  in  the  cen¬ 
tral  Mediterranean  by  more  than  30  Wm-2.  Sensible 
fluxes  (not  shown)  are  comparable  in  magnitude  over 
the  ocean,  with  large  gradients  also  near  the  coast  of 
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(a)  Mean  SW  flux  for  NCEP2  (a)  NCEP2  latent  heat  flux 


(b)  Mean  SW  flux  for  ISCCP 


16  18  20  22  24  26  28  30  32 


Figure  5.  Comparisons  of  shortwave  fluxes  for  2001- 
2004.  (a)  Mean  shortwave  flux  from  (a)  NCEP2  and 
(b)  ISCCP.  Units  are  Wm"2. 


Figure  7.  Time  series  of  latent  heat  flux  for  2001- 
2004  in  the  Central  Mediterranean.  Daily  latent  heat 
flux  from  (a)  NCEP2  and  from  (b)  NCEP2  variables 
used  in  the  COARE  algorithm.  Units  are  Wm~2. 


(c)  Mean  LH  flux  for  NCEP2 


(d)  Mean  LH  flux  for  COARE 
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Figure  6.  Comparisons  of  latent  heat  fluxes  for  2001- 
2004.  (a)  Mean  latent  heat  flux  from  (a)  NCEP2  and 
(b)  COARE.  Units  are  Wm"2. 


Africa.  Interestingly,  the  stronger  mean  shortwave  and 
stronger  mean  latent  heat  fluxes  in  NCEP2  nearly  can¬ 
cel  in  the  net  heat  flux  to  give  net  fluxes  about  the  same 
as  the  sum  of  ISCCP  plus  COARE  fluxes. 

Much  larger  differences  are  seen  when  examining  the 
time  series  of  daily  averaged  latent  heat  flux  at  a  partic¬ 
ular  location,  for  example,  at  34.5°N,  14.5°E  (Figure  7). 
The  standard  deviation  of  latent  heat  flux  is  about  40% 
larger  in  NCEP2  than  in  the  COARE  product.  Maxi¬ 


mum  wintertime  heat  losses  from  the  COARE  flux  are 
about  400 Wm-2,  whereas  NCEP2  maximum  heat  loss 
exceeds  500Wm-2.  Some  reduction  in  the  magnitudes 
could  be  the  result  of  nonlinearity  in  the  bulk  algorithm, 
as  daily  averages  of  the  NCEP2  variables  are  used;  how¬ 
ever,  Jiang  et  al.  [2005]  showed  that  the  nonlinearities 
are  small,  at  least  for  the  tropical  Pacific.  These  differ¬ 
ences  would  be  expected  to  produce  substantial  differ¬ 
ences  in  wintertime  SST.  The  larger  NCEP2  anomalies 
occur  throughout  the  study  region. 

Over  most  of  the  ocean  few  objective  measures  ex¬ 
ist  to  evaluate  absolute  flux  product  accuracy.  How¬ 
ever,  in  the  Mediterranean  Sea  a  recent  study  by  Krah- 
mann  et  al.  [2000]  uses  the  ocean  heat  budget  to  infer 
the  seasonal  cycle  of  net  surface  heat  fluxes  (with  error 
bars)  and  evaluate  several  products.  Temporal  changes 
in  ocean  heat  content  are  the  sum  of  the  net  surface 
heat  fluxes  and  the  ocean  heat  transport  divergence; 
therefore,  to  obtain  an  estimate  of  absolute  net  surface 
heat  fluxes  they  subtract  from  the  time  derivative  of 
the  seasonal  cycle  of  ocean  heat  content  (down  to  the 
bottom  of  the  ocean)  an  estimate  of  the  seasonal  cycle 
of  the  heat  transport  through  the  Strait  of  Gibraltar. 
To  compare  the  heat  flux  products  with  the  estimates 
of  Krahmann  et  al.  [2000]  each  version  of  the  net  fluxes 
was  averaged  over  the  entire  Mediterranean  Sea  and 
from  each  was  extracted  a  single  annual  harmonic  for 
the  four-year  series  (Figure  8).  This  estimate  of  the 
NCEP2  seasonal  net  heat  flux  agrees  well  with  estimates 
from  Krahmann  et  al.  [2000]:  a  maximum  flux  of  166 


Comparison  of  Modeled  SST  in  Mediterranean 


7 


Seasonal  cyde  of  net  surface  heat  flux 


Figure  8.  Annual  cycle  of  net  surface  heat  flux  for 
2001-2004  averaged  over  the  entire  Mediterranean  Sea. 
Net  surface  heat  flux  from  NCEP2  (solid)  and  ISCCP 
plus  COARE  (dash).  Units  are  Wm-2. 

Wm'2,  compared  with  a  range  of  acceptable  values  of 
145  to  170Wm“2,  and  a  minimum  of  about  -167Wm“2, 
compared  with  -160  to  -185Wm-2.  The  combined  IS¬ 
CCP  plus  COARE  seasonal  net  heat  flux  has  a  smaller 
seasonal  range  of  fluxes  than  does  NCEP2.  Although 
the  ISCCP/COARE  maximum  value  of  144 Wm-2  lies 
nearly  within  the  specified  range  of  145  to  170Wm-2,  its 
minimum  value  of  -135Wm-2  is  well  outside  the  Kroh- 
mann  et  al.  [2000]  minimum  range. 

3.  Methodology 

3.1.  Mixed  Layer  Model 

The  model  used  in  this  study  is  the  Price- Weller- 
Pinkel  (PWP)  mixed  layer  model  [Price  et  al .,  1986]. 
The  model  was  run  as  a  column  with  2-m  vertical  res¬ 
olution  to  a  depth  of  about  450  m  at  each  point  in  a 
1°  x  1°  grid,  that  is,  neglecting  input  from  adjacent  grid 
points.  However,  some  experiments  were  performed  to 
estimate  the  effect  of  advection  using  observed  currents 
and  SST  gradients.  The  PWP  model  yields  Ekman  ve¬ 
locities  and  these  are  used  in  some  of  the  experiments 
testing  the  contribution  of  horizontal  advection.  The 
vertical  profile  of  shortwave  irradiance  was  based  on 
values  for  Type  1A  water  (R=0.62  and  7=20).  The 
PWP  model  uses  both  a  gradient  and  bulk  Richardson 
number  for  convection;  the  critical  values  for  overturn¬ 
ing  were  set  at  0.65  and  0.25,  respectively.  Vertical 
diffusivity  was  set  at  5xl0-5  ms-1. 

The  model  was  initialized  using  temperature  and 
salinity  profiles  derived  from  the  World  Ocean  Database 
(2001)  [Conkwright  et  a/.,  2002]  climatology  in  mid- 


October  of  each  year  and  was  then  run  using  daily  forc¬ 
ing  fields  interpolated  to  match  the  6-hr  time  step.  The 
model  was  run  for  periods  of  three  months  beginning  in 
mid-October  of  the  years  2001,  2002,  2003,  and  2004. 
The  fall  was  selected  as  a  time  when  vertical  mixing 
processes  would  have  a  large  impact,  owing  to  strong 
stratification  and  relatively  shallow  mixed  layers.  The 
mid-October  start  time  was  selected  to  avoid  the  shal¬ 
lowest  mixed  layers  in  which  small  errors  in  depth  have 
a  large  effect  on  mixed  layer  temperature. 

Trends  toward  higher  temperatures  and  salinities  in 
the  Mediterranean  Sea  [Cazenave  et  al .,  2001;  Krah- 
mann  and  Schott ,  1998]  made  it  necessary  to  adjust  cli¬ 
matological  values  using  observed  SST.  To  distinguish 
changes  in  the  water  column  from  changes  confined  to 
the  surface,  SST  was  compared  with  climatological  SST 
at  each  grid  point  for  March,  the  month  correspond¬ 
ing  to  the  deepest  mixed  layers.  Observed  tempera¬ 
tures  were  typically  higher  than  climatological  SST  by 
0.5-1. 0°C.  A  single  adjustment  for  the  4-yr  period  was 
made  by  adding  the  mean  March  SST  difference  to  the 
surface  temperature  and  tapering  it  linearly  to  zero  at 
600-m  depth.  A  similar  adjustment  was  made  to  clima¬ 
tological  salinity,  using  a  regression  between  tempera¬ 
ture  and  salinity  to  infer  the  T-S  relationship;  salinity 
adjustments  were  negligible  compared  with  the  seasonal 
cycle  of  salinity  variations.  In  addition  to  the  adjust¬ 
ment  to  climatological  T  and  S  for  a  warming  ocean, 
for  each  model  run  at  each  grid  point,  SST  was  used  to 
adjust  the  initial  temperature  profile,  but  only  down  to 
the  depth  of  the  wintertime  mixed  layer.  Salinity  was 
not  adjusted  further. 

Daily  fields  of  shortwave  flux,  longwave  plus  turbu¬ 
lent  fluxes,  wind  stress,  wind  stress  curl  and  freshwater 
flux  were  used  to  force  the  PWP  model,  after  interpo¬ 
lation  to  the  6-hr  model  time  step.  Two  versions  of 
the  net  surface  heat  flux  were  used:  the  sum  of  all  four 
NCEP2  fluxes  and  the  NCEP2  radiative  fluxes  com¬ 
bined  with  the  NCEP2/COARE  turbulent  fluxes. 

In  the  initial  runs  with  the  prescribed  forcing,  the 
most  noticeable  discrepancies  with  climatological  val¬ 
ues  were  a  tendency  for  the  water  column  to  be  too 
saline  and  the  mixed  layer  to  be  too  deep.  Climatologi¬ 
cal  salinity  shows  a  tendency  to  decrease  (freshening  of 
the  water  column)  presumably  in  response  to  seasonal 
rainfall  over  the  October- January  period,  whereas  the 
NCEP  freshwater  fluxes  are  typically  negative  over  this 
period  (an  excess  of  evaporation).  To  create  a  more  re¬ 
alistic  model  response  to  forcing,  a  simple  bias  correc¬ 
tion  was  made  to  both  the  freshwater  and  heat  fluxes. 
The  bias  corrections  help  maintain  both  a  mixed  layer 
depth  (MLD)  that  more  closely  resembles  climatology 
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and  gives  better  agreement  with  the  observed  seasonal 
trends  of  SST. 

3.2.  Heat  Flux  Corrections 

Two  methods  for  adjusting  the  net  surface  heat  flux 
were  examined  using  the  NCEP2  fluxes.  The  first  ad¬ 
justment  was  based  on  the  difference  between  observed 
and  modeled  SST:  a  straight  line  was  fit  to  the  time  se¬ 
ries  of  the  observed  SST  minus  model  SST.  A  heat  flux 
adjustment  was  estimated  as  qadj  =  mpcp  <  h  >,  where 
m  is  the  tendency  of  the  SST  error  in  °C  s“l,  cp  is  the 
specific  heat  of  seawater,  and  <  h  >  is  the  mean  clima¬ 
tological  MLD.  The  second  method  compared  modeled 
and  climatological  heat  content  (the  vertical  integral  of 
temperature),  to  eliminate  the  effect  of  vertical  entrain¬ 
ment  on  SST.  In  this  case  a  heat  flux  adjustment  was 
estimated  as  qadj  =  ampcp,  where  m  is  now  the  ten¬ 
dency  of  the  heat  content  difference  in  m°C  s_1  and 
a  was  set  to  0.3,  to  allow  for  interannual  departures 
from  climatology.  The  magnitude  of  the  correction  was 
generally  less  than  30Wm“2  for  either  method.  Correc¬ 
tions  based  on  SST  tended  to  be  negative  (less  heat  into 
the  ocean),  whereas  corrections  based  on  heat  content 
tended  to  be  positive.  The  need  to  provide  more  cool¬ 
ing  in  the  mixed  layer  to  match  observed  SST,  while 
the  heat  content  in  the  water  column  was  apparently 
decreasing  too  fast  suggests  that  the  model  mixed  lay¬ 
ers  are  slightly  too  deep. 

Corrections  to  the  surface  heat  fluxes  vary  consider¬ 
ably  (Figure  9)  regionally  and  temporally;  corrections 
shown  were  derived  from  the  SST-based  adjustment. 
The  value  of  the  correction  is  relatively  small  compared 
with  the  net  flux;  corrections  are  typically  less  than 
30Wm-2.  For  the  SST-based  method,  the  flux  correc¬ 
tions  using  both  the  RTG  and  TMI  versions  of  SST 
were  compared  (not  shown);  the  corrections  are  qualita¬ 
tively  similar  and  are  about  the  same  magnitude.  Again 
using  the  SST-based  method,  NCEP2  net  fluxes  were 
compared  with  NCEP2  radiative  plus  NCEP2/COARE 
turbulent  fluxes.  Although  bias  corrections  differed  be¬ 
tween  products  by  as  much  as  30Wm“2,  the  bias  mag¬ 
nitudes  were  again  less  than  about  30Wm“2. 

An  analysis  of  the  heat  budget  in  the  Aegean  Sea 
[Poulos  et  a/.,  1997]  gives  a  climatological  net  surface 
heat  flux  of  approximately  -26Wm~2,  a  small  heat  loss 
to  the  atmosphere  that  is  compensated  by  warm  water 
advection.  In  the  PWP  runs  (in  which  advection  is  ne¬ 
glected)  one  would  expect  a  surface  heat  flux  correction 
of  about  this  size.  However,  the  flux  corrections  in  the 
Aegean  Sea  using  SST  have  no  consistent  sign,  whereas 
the  flux  corrections  based  on  heat  content  are  positive, 
suggesting  that  the  flux  biases  are  not  the  result  of  a 


Figure  9.  Heat  flux  corrections.  Net  surface  heat  flux 
at  each  grid  point  was  adjusted  to  match  the  mean  tem¬ 
perature  tendency  to  observed  values  from  the  RTG 
SST  by  adding  a  constant  offset  to  the  fluxes  from 
NCEP2  for  each  year.  Units  are  Wm"2.  Corrections 
were  typically  less  than  30  Wm-2. 


neglect  of  advection. 

3.3.  Freshwater  Flux  Corrections 

In  contrast  to  the  relatively  small  heat  flux  correc¬ 
tions,  freshwater  fluxes  required  corrections  as  large 
as  the  flux  itself.  On  average  freshwater  fluxes  from 
NCEP2  are  negative  over  the  three-month  model  pe¬ 
riod,  causing  the  water  column  salinity  to  increase, 
whereas  climatological  salinity  decreases  in  the  fall.  An 
analysis  of  climatological  fluxes  in  the  Aegean  Sea  shows 
that  precipitation  exceeds  evaporation  from  November 
through  January  [Poulos  et  ai ,  1997]. 

To  compute  the  freshwater  flux  correction  after  the 
initial  run  with  prescribed  fields,  the  model  salinity  pro¬ 
file  S(z)  at  the  end  of  the  run  is  compared  with  the 
January  climatological  profile,  both  vertically  averaged 
to  the  maximum  climatological  mixed  layer  depth  for 
October- January.  An  empirical  adjustment  is  added  to 
the  mean  freshwater  flux  to  improve  the  match  with 
climatological  salinity.  The  freshwater  flux  corrections 
(Figure  10)  are  generally  positive  (except  in  the  Adri¬ 
atic  Sea),  and  typically  reverse  the  sign  of  the  average 
flux  over  the  three- month  period  from  negative  to  posi¬ 
tive  to  reduce  the  modeled  salinity,  consistent  with  the 
climatological  tendency.  Biases  in  the  salinity  tend  to 
make  the  mixed  layer  too  deep  and  therefore  suppress 
higher-frequency  changes  in  the  surface  temperature, 
particularly  early  in  the  model  run;  these  adjustments 
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Figure  10.  Freshwater  flux  corrections.  Salinity  at 
each  grid  point  was  relaxed  toward  climatological  salin¬ 
ity  by  adding  a  constant  offset  to  the  freshwater  fluxes 
from  NCEP2  for  each  year.  Units  are  my-1. 


generally  make  the  mixed  layer  more  shallow,  consistent 
with  climatological  MLD. 

The  freshening  of  the  water  column  and  the  need 
to  apply  a  freshwater  flux  correction  may  be  owing  in 
part  to  neglect  of  the  inflow  from  adjacent  bodies  of 
water,  as  well  as  from  numerous  rivers.  For  example,  in 
the  Aegean  Sea  it  has  been  estimated  that  evaporation 
would  exceed  precipitation  plus  river  inflow,  on  average, 
but  an  inflow  of  fresh  water  from  the  Black  Sea  creates 
a  net  freshening  (positive  flux)  [Poulos  et  a/.,  1997]. 

After  estimating  both  the  freshwater  and  heat  flux 
corrections,  the  model  was  re-run  with  the  adjusted 
fluxes  (Figure  11).  The  constant  heat  flux  correction 
generally  had  a  negligible  direct  effect  on  the  temper¬ 
ature  tendency,  which  was  dominated  by  variability  on 
time  scales  of  1-2  weeks.  The  freshwater  flux  correction 
consistently  gave  a  shallower  mixed  layer  (closer  to  cli¬ 
matology)  and  improved  the  comparisons  with  observed 
temperature  tendency,  as  discussed  below. 

4.  Comparison  of  Model  and 
Observations 

4.1.  Correlations  with  Observed  Tendency 

To  quantify  model  accuracy,  at  each  grid  point  for 
each  year,  correlations  were  computed  between  temper¬ 
ature  tendency  dT/dt  for  the  model  and  the  observed 
SST  (Figure  12).  Using  the  RTG  SST  product,  which 
covers  more  of  the  study  region,  correlations  are  above 
the  95%  significance  level  (typically  0.35)  over  much  of 


(a)  Model  MLT  (bold).  SST  (thin),  buoy  (dash),  dim  *  (b)  Model  MLD  (bold),  dim  * 


Figure  11.  Observed  and  modeled  variables  at  the 
Santorini  buoy  at  36.3°N,  25.5°E  in  the  fall  of  2002. 
(a)  RTG  SST  (thin  line),  3m  buoy  temperature  (dash), 
and  model  (thick)  and  climatological  mixed  layer  tem¬ 
perature  (*).  (b)  Modeled  (thick)  and  climatological 
mixed  layer  depth  (*).  (c)  As  in  (a),  but  temperature 
tendency  dT/dt. 


Figure  12.  Evaluation  of  model  mixed  layer  tem¬ 
perature.  Correlations  between  temperature  tendency 
dT /dt  from  model  and  from  RTG  SST  for  each  year. 


the  region  in  all  years.  The  region  with  the  lowest  corre¬ 
lations  (mixed  layer  temperature  tendency  least  resem¬ 
bles  that  from  the  RTG  SST)  is  the  central  Mediter¬ 
ranean,  west  of  about  20°E. 

Correlations  vary  from  year  to  year,  suggesting  that 
the  model’s  physics  (or  its  forcing  fields)  may  vary  in 
accuracy.  For  example,  in  the  eastern  Mediterranean, 
correlations  are  above  0.6  in  each  year  except  for  2002. 
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Figure  13.  As  in  Figure  12,  except  TMI  SST. 


A  visual  inspection  of  the  heat  fluxes  and  winds  showed 
no  obvious  anomalies  in  2002.  Correlations  between 
model  and  SST  products  using  the  NCEP2  versus  the 
NCEP2/COARE  net  surface  heat  flux  products  are 
nearly  identical,  suggesting  that  different  flux  products 
are  not  an  important  factor  in  the  differences  in  corre¬ 
lations. 

To  examine  the  possibility  that  the  model  may  be 
performing  well,  but  that  the  SST  used  for  the  evalu¬ 
ation  is  not  uniformly  accurate,  correlations  were  also 
performed  using  the  TMI  SST  in  the  regions  in  which  it 
is  available  (Figure  13).  Again,  correlations  are  signif¬ 
icant  over  much  of  the  region  for  most  years;  however, 
the  regions  and  years  of  high  correlations  differ  from 
those  using  RTG.  Typically,  correlations  of  the  model 
dT /dt  are  higher  with  TMI  than  with  RTG  SST. 

The  RTG  product  is  derived  from  infrared  data, 
which  is  readily  contaminated  by  clouds,  whereas  the 
coarser  microwave  data  has  significant  errors  only  when 
it  is  raining  or  near  land  (as  discussed  above).  A  sim¬ 
ple  indicator  of  the  level  of  cloudiness  can  be  obtained 
from  the  shortwave  radiation  from  ISCCP,  which  de¬ 
creases  substantially  in  the  presence  of  clouds.  The 
effect  of  clouds  on  the  model/SST  comparisons  is  ex¬ 
amined  more  closely  at  one  location  (36.5°N,  18.5°E) 
in  the  central  Mediterranean  (Figure  14).  Correlations 
of  model  dT/dt  with  the  TMI  dT/dt  are  above  0.5  in 
all  years  except  2003,  whereas  the  correlations  with  the 
RTG  dT/dt  are  only  significant  in  2001.  For  this  loca¬ 
tion  the  temperature  tendency  from  each  SST  product 
and  from  the  mixed  layer  model  axe  plotted  in  Figure  14 
(left  column)  for  each  year.  In  the  right  column  are  plot¬ 
ted  the  shortwave  flux  anomalies  from  the  seasonal  cy¬ 
cle.  There  are  relatively  large  negative  anomalies  (large 


Effect  of  clouds  on  RTG  SST 
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Figure  15.  Effect  of  cloudiness  on  temperature  ten¬ 
dency.  Scatter  plot  of  correlations  between  tempera¬ 
ture  tendency  from  RTG  and  from  microwave  SST  (as  a 
proxy  for  RTG  accuracy)  versus  shortwave  flux  anoma¬ 
lies  (as  a  proxy  for  cloud  cover).  The  regression  line 
(dashed)  between  the  proxy  for  RTG  SST  accuracy  and 
the  proxy  for  cloud  cover. 


cloud  contributions)  except  in  year  2001,  the  only  year 
in  which  RTG  compares  well  with  the  model.  This  ex¬ 
ample  suggests  that  cloudiness  may  be  responsible  for 
lower  correlations  between  the  model  and  RTG  dT/dt. 

The  effect  of  cloudiness  on  infrared  SST  tendency  can 
be  estimated.  If  the  clouds  are  significantly  degrading 
the  RTG  SST  (but  not  the  microwave  SST),  one  would 
expect  poor  correlations  between  the  two  SST  prod¬ 
ucts  during  cloudy  periods  (Figure  15).  For  the  cen¬ 
tral  Mediterranean  (33.5-37. 5°N,  15.5-21.5°E)  where 
the  apparent  model  accuracy  differs  greatly  between 
the  two  SST  products,  correlations  between  dT /dt  from 
RTG  and  TMI  are  plotted  against  the  average  anoma¬ 
lies  of  shortwave  flux  for  15  October  -  15  December 
of  each  year.  From  the  plot  it  can  be  seen  that  low 
correlations  of  RTG  with  TMI  SST  coincide  with  neg¬ 
ative  flux  anomalies  (cloudiness),  indicating  that  errors 
from  cloud  contamination  in  the  RTG  product  are  sig¬ 
nificantly  reducing  its  usefulness  for  evaluating  model 
performance. 

Thus,  the  TMI  SST  is  a  better  indicator  of  model 
performance  than  RTG,  where  the  data  are  available 
and  are  not  influenced  by  nearby  land.  In  the  Aegean 
Sea,  where  TMI  data  are  masked  out  by  land,  the  model 
performs  quite  well  in  all  four  years,  suggesting  that 
cloud  contamination  of  RTG  SST  may  be  less  of  a  prob- 
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(b)  dT/dt  2002 
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Figure  14.  Temperature  tendency  and  shortwave  flux  anomalies  at  36.5°N,  18.5°E.  Temperature  tendency  from 
RTG  (thin  line),  from  microwave  SST  (thick  line),  and  from  the  mixed  layer  model  (dashed  line)  for  the  fall  of  (a) 
2001,  (b)  2002,  (c)  2003,  and  (d)  2004.  (e)-(h)  Anomalies  from  seasonal  cycle  of  shortwave  radiation  for  fall  of  the 
same  four  years. 
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lem  there. 

4.2.  Mixed  Layer  Temperature  Variability 

The  PWP  model  shows  substantial  variations  in 
mixed  layer  temperature  (MLT)  variability  on  time 
scales  of  about  10-14  days.  (Recall  that  variability  on 
time  scales  shorter  than  about  9  days  has  been  removed 
from  the  SST  data,  as  described  in  Section  2.3).  These 
changes  could  be  caused  by  variations  in  heat  flux,  in 
wind  forcing,  or  mixing.  In  addition,  the  column  mixed 
layer  model  neglects  temperature  advection  from  both 
Ekman  and  geostrophic  current  components. 

To  determine  the  extent  of  the  various  processes  to 
dT/dt ,  each  of  the  forcing  functions  in  turn  was  low- 
pass  filtered  (half- power  of  30  days),  while  the  others 
retained  daily  variations.  A  location  and  year  in  which 
the  correlations  with  both  RTG  and  TMI  were  rela¬ 
tively  high  was  selected:  33.5°N,  27.5°E  in  2003.  The 
correlation  between  model  and  observations  was  0.65 
for  the  baseline  run  (Figure  16a).  Removing  the  high- 
frequency  fluctuations  in  wind  stress,  wind  stress  curl, 
or  shortwave  radiation  had  a  negligible  effect  on  the  cor¬ 
relations  with  TMI.  However,  removing  high-frequency 
fluctuations  from  the  turbulent  fluxes  (actually  the  tur¬ 
bulent  plus  longwave  fluxes,  but  longwave  fluxes  are 
quite  small)  removed  nearly  all  the  variability  in  the 
model  dT/dt ,  reducing  the  correlation  to  0.39  (Fig¬ 
ure  16b). 

The  comparison  of  the  model  and  observed  dT/dt  in 
Figure  16a  is  fairly  typical  in  that  the  model  slightly 
underestimates  the  magnitudes,  even  at  the  beginning 
of  the  run  when  model  MLD  agrees  well  with  climatol¬ 
ogy.  Near  the  end  of  the  run,  at  many  locations,  MLD 
is  overestimated,  resulting  in  even  smaller  dT/dt  mag¬ 
nitudes,  relative  to  observed  values.  These  comparisons 
suggest  that  the  energetic  NCEP2  turbulent  fluxes  are 
more  consistent  with  observed  SST  changes  than  the 
NCEP2/COARE  fluxes. 

The  PWP  model  run  independently  at  each  grid 
point  neglects  horizontal  processes  (advection  and  dif¬ 
fusion)  that  might  significantly  affect  the  mixed  layer 
temperatures.  Experiments  were  performed  to  esti¬ 
mate  the  effects  of  temperature  advection  in  the  mixed 
layer.  To  estimate  the  effect  of  horizontal  advection 
by  the  Ekman  response  to  winds,  the  modeled  Ekman 
velocity  was  combined  with  observed  SST  gradients  (as¬ 
suming  a  uniform  horizontal  temperature  gradient  with 
depth  over  the  Ekman  layer)  and  compared  with  the 
model  run  without  advection.  The  difference  at  sev¬ 
eral  grid  points  showed  a  relatively  small  contribution, 
a  warming  or  cooling  of  about  0.1-0.2°C  over  the  3- 
month  period.  Correlations  between  modeled  and  ob- 


(a)  Unfiltered  turbulent  flux  (dashed) 


Figure  16.  Source  of  high  frequency  variations  in 
mixed  layer  temperature.  Temperature  tendency  dT/dt 
at  33.5°N,  27.5°E  in  fall  of  2003  from  (a)  TMI  SST 
(solid  line)  and  model  run  with  daily  forcing  (dashed 
line)  and  from  (b)  model  run  with  smoothed  version 
of  turbulent  heat  fluxes  (dashed  line).  Observed  SST 
repeated  in  (b)  (solid  line).  Units  are  °Cs_1. 


served  dT/dt  were  not  significantly  changed.  The  PWP 
model  does  not  include  geostrophic  currents.  Hypothet¬ 
ically,  these  could  be  estimated  using  geostrophic  veloc¬ 
ity  anomalies  from  the  altimeter;  however,  a  mean  sea 
surface  for  the  Mediterranean  Sea  is  not  readily  avail¬ 
able  to  supplement  the  anomalies. 

Currents  were  available  at  3-m  depths  at  the  buoys 
in  the  Aegean  Sea,  although  there  were  substantial  data 
dropouts  during  the  study  period.  A  nearly  complete 
record  for  the  fall  of  2002  at  the  Santorini  buoy  (Fig¬ 
ure  17)  was  used  to  examine  the  effect  of  Ekman  plus 
geostrophic  advection.  The  observed  daily-averaged 
currents  resemble  the  modeled  Ekman  currents,  but  the 
time  series  are  only  marginally  correlated.  The  mod¬ 
eled  currents  do  not  have  the  observed  lower  frequency 
(>  10-day  periods)  fluctuations,  which  presumably  are 
the  geostrophic  component.  The  effect  of  advection  on 
model  accuracy  by  the  observed  currents  is  quantified 
(Table  2)  using  correlations  between  model  and  3-m 
buoy  temperature  tendency  and  between  the  model  and 
RTG  SST  tendency.  Adding  advection  increased  (buoy) 
or  decreased  (RTG)  correlations  somewhat,  but  these 
changes  are  relatively  small.  Again,  this  shows  that  the 
PWP  column  model  is  doing  fairly  well  in  hindcasting 
MLT. 
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(a)  Zonal  velocity 


(a)  Std  dev  of  QuikSCAT  wind  speed 
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(b)  Meridional  velocity 


Figure  IT.  Observed  and  modeled  currents  at  San¬ 
torini. (a)  zonal  and  (b)  meridional  currents  from  buoy 
(dashed)  and  model  (solid). 


Figure  18.  Standard  deviation  of  wind  speeds  for 
November-December-January.  (a)  QuikSCAT  and  (b) 
NCEP2  in  units  of  ms-1. 


(b)  Std  dev  of  NCEP2  wind  speed 


Table  1.  Correlations  of  dT/dt  with  and  without  ad- 
vection  at  Santorini 


RTG  SST 

Buoy  3-m 

No  advection 

0.59 

0.60 

Advection 

0.51 

0.63 

4.3.  Effect  of  Wind  Accuracy 

As  noted  in  Section  4.2,  wind  stress  and  wind  stress 
curl  have  little  direct  effect  on  the  prediction  of  SST 
using  the  model,  compared  with  the  effects  of  turbulent 
heat  fluxes.  However,  wind  speed  likely  makes  a  large 
contribution  to  the  NCEP2  turbulent  fluxes,  and  poor 
spatial  resolution  in  wind  speed  may  have  a  large  ef¬ 
fect  on  SST  through  the  turbulent  fluxes.  Wind  speeds 
near  the  numerous  islands  in  the  Mediterranean  Sea  are 
highly  variable  spatially,  owing  to  the  complex  topog¬ 
raphy,  and  are  not  resolved  by  the  NCEP2  Reanalysis 
product. 

To  study  the  effect  of  poor  spatial  resolution  in  wind 
speed  on  the  NCEP2  turbulent  flux  products,  and  in 
turn  on  model  accuracy,  a  region  south  of  the  island 
of  Crete  was  selected  for  further  study.  High-resolution 
fields  of  wind  speed  derived  from  the  QuikSCAT  data 
(Figure  lb)  show  variability  on  spatial  scales  much 
smaller  than  the  resolution  of  the  NCEP2  flux  prod¬ 
ucts;  these  differences  are  sufficiently  robust  that  they 
can  be  seen  in  the  standard  deviation  of  the  November- 
December-January  wind  speeds  for  2001-2004  (Fig¬ 
ure  18).  Note  the  local  minimum  in  wind  speed  vari¬ 
ations  south  of  central  Crete  (35°N,  25°E)  where  the 


topography  blocks  the  prevailing  northerly  winds.  The 
QuikSCAT  wind  fields  used  here,  although  mapped  to 
the  same  1°  grid  as  the  NCEP2  products,  have  true  1° 
resolution,  whereas  the  winds  (and  the  other  variables) 
used  in  the  NCEP2  and  NCEP/COARE  products  are 
much  smoother. 

To  characterize  the  accuracy  of  the  NCEP2  fluxes, 
correlations  were  performed  between  the  NCEP2  and 
QuikSCAT  winds  along  34.5°N  from  19.5-32.5°E  for 
each  of  the  four  years  of  the  study,  on  the  assump¬ 
tion  that  locations  where  NCEP2  winds  are  more  highly 
correlated  with  the  QuikSCAT  winds  have  more  accu¬ 
rate  flux  estimates  than  locations  with  low  correlations. 
This  proxy  for  flux  field  accuracy  was  then  compared 
with  a  measure  of  model  accuracy,  the  correlation  be¬ 
tween  observed  and  model  temperature  tendency  dT /dt 
(Figure  19).  A  regression  between  the  two  correlations 
gives  the  best-fit  line  (dashed  line  in  Figure  19)  with  a 
correlation  of  0.41,  which  is  just  significant  at  95%  con¬ 
fidence.  This  analysis  gives  an  estimate  of  the  degra¬ 
dation  in  model  performance  (a  drop  in  correlations  as 
large  as  0.2)  that  results  from  poorly  resolved  spatial 
variations  in  the  flux  products. 

4.4.  Characterizing  Model  Accuracy 

In  Section  4.1  the  correspondence  between  model  and 
observed  temperature  tendency  was  characterized  using 
correlations.  One  drawback  of  the  use  of  a  correlation 
is  that  it  does  not  include  any  measure  of  the  relative 
magnitudes  of  the  model  and  the  observed  fluctuations, 
and  therefore,  hides  an  under-  or  over-prediction.  For 
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Effect  of  wind  speed  accuracy  on  model  accuracy 
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Correlation  of  NCEP2  and  QuikSCAT  wind  speed 

Figure  19.  Effect  of  wind  speed  accuracy  on  model 
accuracy.  Correlations  between  observed  and  mod¬ 
eled  temperature  tendency  dT/dt  versus  correlations 
between  NCEP2  and  QuikSCAT  wind  speed  south  of 
Crete. 

example,  the  difference  in  the  magnitudes  of  the  latent 
heat  flux  anomalies  of  NCEP2  versus  NCEP2/COARE, 
as  shown  in  Figure  7,  along  with  the  sensitivity  of  the 
model  to  latent  heat  flux  anomalies,  suggests  that  the 
magnitude  of  the  model’s  response  should  be  larger  for 
NCEP2.  Correlations  with  observed  temperature  ten¬ 
dency  for  the  model  forced  with  the  two  flux  products 
were  nearly  identical  and  give  no  indication  of  which 
product  gives  smaller  errors.  An  additional  problem 
is  that  a  correlation  between  model  and  observations 
can  be  increased  by  filtering  the  model  input  (or  out¬ 
put)  to  remove  higher-frequency  (uncorrelated)  vari¬ 
ations  without  actually  improving  the  model  perfor¬ 
mance.  The  filtering  may  result  in  under-prediction  by 
the  model,  but  this  effect  will  not  be  measured  with  a 
correlation.  Thus,  a  correlation  is  not  an  ideal  method 
of  evaluation. 

An  alternative  method  to  describe  model  accuracy  is 
the  so-called  Taylor  diagram  [ Taylor ,  2001].  Polar  coor¬ 
dinates  are  used  to  plot  both  correlation  (0)  and  magni¬ 
tudes  (r).  Here,  the  normalized  version  of  this  diagram 
is  used,  so  that  the  root-mean-square  (rms)  ratio  of  the 
modeled  values  to  the  observed  values  is  plotted  as  the 
radial  distance  from  the  origin  (Figure  20).  The  angle 
with  respect  to  the  x-axis  is  derived  from  the  correla¬ 
tion,  as 

9  —  cos~l p 

so  that  a  perfect  correlation  lies  along  the  x-axis  and  a 


Figure  20.  Example  of  a  normalized  Taylor  diagram. 
Using  polar  coordinates,  the  ratio  of  the  standard  de¬ 
viations  of  model  and  observed  temperature  tendency 
(radial  coordinate)  and  the  correlation  angle  6  are  plot¬ 
ted  for  each  location  and  each  year.  The  normalized 
error  is  the  distance  r  from  each  point  to  the  location 
of  a  perfect  correlation  (on  the  x-axis)  and  a  ratio  of 
one.  See  text  for  explanation. 


zero  correlation  lies  along  the  y-axis.  The  center  of  the 
cluster  of  points  in  the  figure  lies  at  a  magnitude  ratio 
of  about  0.6  and  correlations  of  about  p  =  0.66. 

The  distance  of  any  point  from  the  intersection  of  the 
circle  of  radius  one  and  the  x-axis  is  a  normalized  mea¬ 
sure  of  the  accuracy  of  the  model  prediction  of  temper¬ 
ature  tendency;  the  distance  is  the  ratio  of  the  standard 
deviation  of  the  error,  divided  by  the  standard  devia¬ 
tion  of  the  signal,  here,  the  observed  value  of  dT /dt.  In 
Figure  20  the  error  r  shown  is  about  0.75,  which  means 
at  this  location  in  this  year,  the  PWP  model  has  an  er¬ 
ror  of  about  75%  of  the  standard  deviation  of  observed 
dT/dt.  For  r  >  1  the  model  error  would  exceed  the 
signal  and  the  model  has  no  useful  skill  in  simulating 
dT/dt. 

Normalized  errors  using  COARE  versus  NCEP2  heat 
fluxes  (not  shown)  are  approximately  10%  more  in  the 
central  Mediterranean,  a  clear  indication  of  an  under- 
prediction  of  temperature  tendency  by  COARE  there. 
The  use  of  dT /dt  as  a  metric  for  evaluating  the  model  is 
quite  stringent  and,  therefore  the  normalized  distances 
are  large.  Clearly  the  model  simulates  MLT  well  (for 
example,  Figure  I  la)  and  a  metric  based  on  T(t)  would 
likely  have  much  smaller  normalized  errors. 

The  measurements  of  error  in  Figure  21  use  the  SST 
product  that  gives  the  smallest  normalized  errors  over 
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Figure  21.  Evaluation  of  model  using  Taylor  dia¬ 
gram.  Normalized  errors  in  temperature  tendency  by 
the  model  relative  to  the  best  SST  observation. 


all  four  years  (compare  with  Figure  12  and  13).  Gen¬ 
erally,  the  TMI  SST  is  used  for  the  interior  points  and 
RTG  SST  is  used  near  land.  A  distance  of  1.0  indicates 
where  the  expected  model  error  is  the  same  size  as  the 
standard  deviation  of  dT / dt ,  that  is,  where  the  PWP 
model  has  no  useful  skill  in  predicting  MLT  variations. 
The  regions  of  largest  (relative)  errors  are  the  central 
Mediterranean  Sea,  the  western  part  of  the  study  re¬ 
gion. 


5.  Summary  and  Conclusions 

The  goal  of  this  study  is  to  determine  whether  avail¬ 
able  observations  can  be  used  to  determine  when  an 
ocean  model  is  missing  essential  physics,  specifically, 
whether  mixing  processes  not  represented  in  a  mixed 
layer  model  are  responsible  for  systematic  temperature 
prediction  biases.  At  this  point,  errors  in  forcing  fields 
and  comparison  temperature  data  make  that  determi¬ 
nation  difficult.  Based  on  the  analyses  here,  the  model 
appears  to  be  performing  well  compared  with  the  avail¬ 
able  observations,  in  that  times  and  locations  of  poor 
model  performance  generally  correspond  to  poor  qual¬ 
ity  of  forcing  or  observed  SST  fields.  For  example, 
poor  model  performance,  as  judged  by  correlations  with 
temperature  tendency  from  the  RTG  SST  product,  cor¬ 
respond  to  periods  of  heavy  cloud  cover  or  to  regions 
where  spatial  resolution  in  air-sea  fluxes  is  inadequate. 

On  the  other  hand,  it  appears  that  the  model  can 
help  determine  which  of  the  available  observations  and 
forcing  fields  are  most  accurate.  For  example,  the  mi¬ 
crowave  SST  gave  consistently  higher  correlations  with 


the  model  than  the  infrared- based  RTG  product,  ex¬ 
cept  near  land  where  there  is  a  warm  SST  bias  in  the 
microwave  data.  The  accuracy  of  the  infrared  product 
is  degraded  in  periods  of  high  cloud  cover  (Figure  15), 
as  parameterized  by  negative  shortwave  flux  anomalies. 
Large  corrections  to  the  freshwater  fluxes  (in  the  sense 
of  increasing  the  precipitation)  are  required  for  consis¬ 
tency  with  climatological  salinity,  suggesting  that  these 
fluxes  have  little  useful  skill  for  forcing  the  model. 

NCEP2  heat  fluxes,  which  have  the  largest  temporal 
variations  of  the  flux  fields  examined  here,  still  under¬ 
estimate  the  variations  in  temperature  tendency  (Fig¬ 
ure  16).  Part  of  the  underestimate  comes  from  an 
overestimate  of  mixed  layer  depth;  however,  the  un¬ 
derestimate  is  apparent  even  at  the  beginning  of  the 
model  runs,  when  the  MLD  is  very  near  its  climatolog¬ 
ical  value  and  occurs  at  locations  where  MLD  matches 
climatology  well  throughout  the  run,  as  in  Figure  16. 
NCEP2  radiative  fluxes  have  temporal  variations  about 
twice  as  large  those  from  ISCCP  and  temporal  varia¬ 
tions  in  turbulent  fluxes  about  40%  larger  those  made 
by  using  NCEP2  daily  fields  in  the  COARE  algorithm. 
The  COARE  algorithm  is  designed  to  produce  realistic 
fluxes  using  hourly  input  fields,  rather  than  the  daily 
values  used  here,  which  could  account  for  the  under¬ 
estimate;  however,  a  recent  analysis  by  [Jiang  et  al., 
2005]  suggests  that  the  COARE  algorithm  is  not  highly 
nonlinear. 

The  source  of  the  high-frequency  temperature  fluc¬ 
tuations  was  determined  using  a  series  of  model  experi¬ 
ments  and  comparing  the  results  with  observations.  At 
the  highest  frequency  resolved  by  the  observed  SST  (pe¬ 
riods  of  about  9  days)  variations  in  SST  are  caused  by 
corresponding  variations  in  turbulent  heat  fluxes.  No 
other  forcing  field  examined  (shortwave  radiation,  wind 
stress,  freshwater  fluxes)  makes  a  significant  contribu¬ 
tion.  Advection  makes  only  a  small  contribution  to 
the  SST  tendency  in  the  Aegean  Sea,  the  only  location 
where  currents  axe  available. 

As  expected  from  the  importance  of  the  turbulent 
heat  flux,  errors  in  the  flux  estimates  appear  to  de¬ 
grade  model  performance.  In  the  region  south  of  Crete, 
where  QuikSCAT  winds  show  large  topographic  effects 
at  small  scales,  the  coarse  resolution  of  NCEP2  flux 
fields  appears  to  significantly  degrade  model  results. 
Using  the  correlation  of  NCEP2  and  QuikSCAT  wind 
speeds  (Figure  19)  as  a  measure  of  NCEP2  turbulent 
flux  accuracy,  regions  of  poor  NCEP2  flux  accuracy  cor¬ 
respond  to  regions  of  poor  model  performance. 

Model  accuracy  is  characterized  using  correlations  of 
observed  and  modeled  temperature  tendency,  to  em¬ 
phasize  the  highest  frequencies  resolved  by  the  observa- 
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tions  (periods  of  9  days  or  more).  To  include  informa¬ 
tion  about  over-  or  under-prediction  of  the  temperature 
tendency,  normalized  errors  are  computed  based  on  the 
Taylor  diagram.  This  characterization  reveals  that  the 
PWP  model  (forced  by  currently  available  fields)  has 
its  greatest  skill  in  predicting  temperature  tendency  in 
the  eastern  part  of  the  study  region,  and  has  no  use¬ 
ful  skill  in  the  central  Mediterranean  (15°E-  20°E),  the 
westernmost  part  of  the  study  region.  Further,  it  shows 
that  the  larger  NCEP2  flux  anomalies  give  better  re¬ 
sults  than  the  other  flux  products  tested,  particularly 
in  the  eastern  Mediterranean.  Skill  in  modeling  temper¬ 
ature,  rather  than  temperature  tendency,  is  expected  to 
be  much  greater. 

The  new  satellite  products  (microwave  SST  and  scat- 
terometer  winds)  clearly  have  the  potential  to  improve 
both  forcing  and  comparison  fields,  to  allow  an  exami¬ 
nation  of  shortfalls  in  model  physics.  However,  consis¬ 
tent  high-spatial-resolution  SST  and  flux  products  need 
to  be  derived  from  the  new  data  fields  to  allow  inter¬ 
pretation  of  poor  model  performance  as  missing  or  in¬ 
accurate  model  physics.  A  new  SST  product  (GHRSST 
from  GODAE)  combining  infrared,  microwave,  and  in 
situ  temperatures  will  soon  be  available,  as  will  a  higher 
resolution  (12.5km)  wind  product  from  the  QuikSCAT 
scatterometer.  As  suggested  by  the  analysis  here,  these 
products  may  make  possible  the  detection  of  shortcom¬ 
ings  in  model  physics  by  careful  comparison  with  obser¬ 
vations. 
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